* Figure4stderr.do  02/04/2024
* Figure 4 standard errors
* expanded file has the full population of opioid deaths, and a random sample of U.S. residents not dying from opioids

do Figure4prep
keep if ucd=="Drug" & mcd=="Opioid" & year>=2012 & (age=="0-44" | age=="45-64" | age=="65+")
drop if deaths==.
gen byte black = race==`"{"Black or African American"}"'
gen byte female = gender==`"{"Female"}"'
encode age, gen(agenum)
encode state, gen(censusdiv)
local kprt=40                                  // to economize on memory, cut the survivors' dataset by a factor of `kprt'
gen pop1 = deaths
gen pop0 = round((pop - deaths)/`kprt')        // Did not die from opioids during the year
drop deaths pop
reshape groups diedopioid 0 1
reshape cons censusdiv female agenum black year
reshape var pop
reshape long
compress
replace diedopioid = diedopioid*100000/`kprt'  // Deaths per 100K, with `kprt' adjustment for oversample
drop if pop==0                                 // No synthetic individuals from empty cells
expand pop
drop pop

gen gapraw = .
gen gapadj = .
gen gapold = .
gen gapyoung = .
gen gapmid = .
gen gapeld = .
gen seadj = .
gen seold = .
gen seyoung = .
set more off
foreach yr of numlist 2012/2021 {
  disp `yr'
  reg diedopioid black if year==`yr'
  replace gapraw = _b[black] if year==`yr'
  reg diedopioid black female i.agenum i.censusdiv if year==`yr'
  replace gapadj = _b[black] if year==`yr'
  replace seadj = _se[black] if year==`yr'
  reg diedopioid black female i.agenum i.censusdiv if year==`yr' & agen!=1
  replace gapold = _b[black] if year==`yr'
  replace seold = _se[black] if year==`yr'
  reg diedopioid black female i.censusdiv if year==`yr' & agen==1
  replace gapyoung = _b[black] if year==`yr'
  replace seyoung = _se[black] if year==`yr'
  reg diedopioid black female i.censusdiv if year==`yr' & agen==2
  replace gapmid = _b[black] if year==`yr'
  reg diedopioid black female i.censusdiv if year==`yr' & agen==3
  replace gapeld = _b[black] if year==`yr'
}
gen gapadjupper = gapadj + 1.96*seadj
gen gapadjlower = gapadj - 1.96*seadj
gen gapoldupper = gapold + 1.96*seold
gen gapoldlower = gapold - 1.96*seold
gen gapyoungupper = gapyoung + 1.96*seyoung
gen gapyounglower = gapyoung - 1.96*seyoung
gen zeroes = 0
sort year censusdiv female agenum black
gen byte kp = year!=year[_n-1]
keep if kp>0
#delimit ;
twoway (line zeroes year, lcolor(gray) lpattern(dash) xlabel(2012[1]2021))
  (rcap gapadjlower gapadjupper year, lcolor(black))
  (rcap gapoldlower gapoldupper year, lcolor(green))
  (rcap gapyounglower gapyoungupper year, lcolor(blue))
  (connected gapadj gapold gapyoung year, lcolor(black green blue) mcolor(black green blue) msize(small small small)
  title("Figure 4.  Black-white gaps in opioid fatality rates") subtitle("adjusted for gender, age group, and Census division")
  legend(order(6 "Ages 45+" 5 "All ages" 7 "Ages 0-44")  region(lcolor(white)))
  graphregion(color(white)) ylabel(,nogrid) xtitle("") ytitle("Deaths per 100K population")
  note("Notes: For each year an opioid death indicator is regressed on indicators for race, gender, age group, and Census division, with"
  "100K times the black coefficient shown in the figure.  A regression observation is a black or white U.S. resident alive January 1."
  "Population and mortality source: CDC Wonder.", size(vsmall)) );
#delimit cr
label data "Confidence intervals from year-by-year individual-level opioid-mortality regressions"
save Figure4stderr
